Grass demonstration

In this lab, we will be introduced to another open source GIS software named Grass. GRASS (Geographic Resources Analysis Support System) has been under continuous development since 1982 and during the 90s a large number of federal US agencies, universities, and private companies were involved in the development. The core components of GRASS 1.0-5.0 and the management of the integration efforts into GRASS GIS releases were accomplished by the U.S. Army - Construction Engineering Research Laboratory (USA-CERL) in Champaign, Illinois. USA-CERL completed its last release of GRASS as version 4.1 in 1992, and provided five updates and patches to this release through 1995. USA-CERL also wrote the core components of the GRASS 5.0 floating point version. In 1997 and after two years of uncertainty, GRASS development was taken over by academia. Since then the international GRASS Development Team manages the source code, releases and documentation.

One of the advantages of GRASS GIS is its numerous modules for analyzing all manners of spatial and temporal data. GRASS GIS has over 500 different modules in the core distribution and over 300 addons

Data and software

Included in this data are LAS files, reference maps for those files and an orthophoto of Washtenaw county. A LAS file (.las) is the standard binary format for storing airborne lidar data. The LAS dataset allows you to examine LAS files, in their native format, quickly and easily, providing detailed statistics and area coverage of the lidar data contained in the LAS files. The reference map indicates the specific location of the .las in the county. We will drape the orthophoto over our interpolated DSM.

Note: On the Mac OS, the r.in.lidar module may encounter problems. Mac user may need to use r.in.ascii in combination with libLAS instead of the provided code. 3D view may be inaccessible for Macs.

Basic introduction to graphical user interface

GRASS GIS Spatial Database

GRASS uses specific database terminology and structure (GRASS GIS Spatial Database):

  • A GRASS GIS Spatial Database (GRASS database) consists of a directory with specific Locations (projects) where data (data layers/maps) are stored.

  • Location is a directory with data related to one geographic location or a project. It is important to remember and this is a benefit of GRASS that all data within one Location has the same coordinate reference system.

  • A Mapset is a collection of maps within a Location, containing data related to a specific task, user or a smaller project.

Exercise instructions

  1. Start GRASS GIS. Now, you will create a new directory where you data will be stored using the interface - hit New. Open the Location Wizard with the New button on the left part of the welcome screen. Select a name for the new location(this is where all your data will be stored)and create this folder.

alt text alt text

We will establish this location’s projection-CRS (coordinate reference system). Our LiDAR data uses NAD83(HARN) / Michigan South. We will search for the correct EPSG using the Select EPSG code of spatial reference system. Use the search tool to find the correct EPSG(). After choosing a projection, a new location will be listed on the start-up screen.

alt textalt text

Now, select the new location and mapset PERMANENT and press Start GRASS session.

Importing data

In this step we will import data. This is an important step in GRASS, as it uses a specific file format only readibly in the software. Using the menu choose File - Import raster data select Common formats import In the resulting dialog box, browse to find the DEM file, change the name to DEM, and click the button Import. The imported layers will be added to GUI automatically (you can add them manually using alt text and alt text.

You can also use the command line:r.import input=???.tif output=

Computational region

It important to understand the computational region in Grass, which confines any analysis to a specific geographic area. The computational region can be set, subsetting a larger data extent, for quicker testing or analysis of specific regions based on administrative units. You must always set the computational region properly. These are some things keep in mind when using the computational region function:

  • it is defined by region extent and raster resolution
  • it applies to all raster operations
  • it is persistent between GRASS sessions, but can be different for different mapsets
  • advantages: keeps your results consistent, avoids clipping, for computationally demanding tasks set region to smaller extent, offers a simple way to check your result before running your analysis on a larger computational region (computationally expensive operation).

alt text alt text alt text

You can check the computational region in the Console using:

g.region -p

The computational region can be set using a raster map:

g.region raster=Arb_DEM -p

Resolution can be set separately using the res parameter of the g.region module. The units are the units of the current location, in our case meters. This can be done in the Resolution tab of the g.region dialog or in the command line in the following way (using also the -p flag to print the new values):

g.region res=3 -p

The new resolution may be slightly modified in this case to fit into the extent which we are not changing. However, often we want the resolution to be the exact value we provide and we are fine with a slight modification of the extent. That’s what -a flag is for. On the other hand, if alignment with cells of a specific raster is important, align parameter can be used to request the alignment to that raster (regardless of the extent).

The following example command will use the extent from the raster named ortho, use resolution 5 meters, modify the extent to align it to this 5 meter resolution, and print the values of this new computational region settings:

g.region raster=ortho res=5 -a -p

Modules

Different functions are available through modules in GRASS (which are sometimes called tools, functions, or commands). Modules respect the following naming conventions:

Prefix Function Example
r. raster processing r.mapcalc: map algebra
v. vector processing v.surf.rst: surface interpolation
i. imagery processing i.segment: image segmentation
r3. 3D raster processing r3.stats: 3D raster statistics
t. temporal data processing t.rast.aggregate: temporal aggregation
g. general data management g.remove: removes maps
d. display d.rast: display raster map

Modules staring with v.net. deal with vector network analysis. The name of the module helps to understand its function, for example v.in.lidar starts with v so it is a vector maps function, the following name indicates module functionality importing, and finally lidar indicates that it deals with lidar point clouds. Modules and their descriptions with examples can be found the documentation. The documentation is included in the local installation and is also available online.

Displaying raster and vector data

Here is a good tutorial showing how to manipulate and display your maps link

LiDAR

The fastest way to analyze a LiDAR point cloud is to use binning and create a raster map. In Grass, the r.in.lidar converts a point count (point density) into raster. As we don’t know the spatial extent of the point cloud, we cannot set the computation region ahead. Instead, we can tell the r.in.lidar module to determine the extent automatically using the -e flag. We can set the resolution (pixel size). Remember that the resolution impact the processing speed. The -r flag sets the computational region to match the output. The -o flag ignores the projection definition (the data do not have a projection, but we know that it is NAD83(HARN) / Michigan South)

r.in.lidar input=297285v.las output=Arb_297285 method=n -e -n -o resolution=10

After it has been successfully imported, we can see the spatial patterns of the LiDAR data. From the LiDAR metadata we know that the data is in feet. Grass automatically uses meters and so if you set the resolution to 3, it results in pixel size of 0.9144. We can also examine variation of the raster using the r.info function:

r.info map=Arb_297285

We can do a quick check of some of the values using the query tool in the Map Display. Examine some the pixel values

Since there are numerous points per one cell, we can use a finer resolution (notice the overwrite flag). This will result in a meter resolution map:

r.in.lidar input=297285v.las output=Arb_297285 method=n -e -n -o resolution=3 --overwrite

We can examine the distribution of the values using the histogram function. Histogram is accessible from the context menu of a layer in Layer Manager, from theMap Display`` toolbarAnalyze map``` button; or you can simply use this function in the console:

d.histogram map=Arb_297285

Now we need to do a bit of housekeeping to ensure that subsequent calculations inherit the same extent and resolution. We do this by applying the extent and resolution of our calculated image using the g.region function:

g.region raster=Arb_297285 -p

However, while this region has a lot of cells, some areas are empty due to scattering of data. To account for this absence of data in some locations, we could modify the resolution of the computational region. We use the a flag -ap to align region to resolution. For this exercise we will maintain the 1 meter resolution:

g.region res=3 -ap

DSM

A Digital Surface Model (DSM) represents the elevations of the reflective surfaces of trees, buildings, and other features elevated above the “Bare Earth”. We can use binning to obtain the DSM. Because we set the g.region, our new data will inherit the resolution and extent from the existing computational region:

r.in.lidar input=297285v.las output=binned_dsm method=max -o

To get a better understanding of different cloud values, we can compute statistics using r.report (Report raster statistics in the layer context menu):

r.report map=binned_dsm units=c

This table indicate that there are a few really high numbers (outliers at around 970 m). To better visualize this we change the color classes based on equalizing the histogram (-e flag). This allows us to see the contrast in the rest of the map (using the viridis color table):

r.colors map=binned_dsm color=elevation -e

We can also identify possible low outliers using the minimum function (method=min, we already used method=max for DSM):

r.in.lidar input=297285v.las output=minimum -o method=min

If we compare this with the range and mean of the low values, it will give an idea of how best to confine the range to eliminate possible outliers.

r.report map=minimum units=c

Again, to see the data, we can use histogram to equalized color table:

r.colors map=minimum color=elevation -e

The zrange parameter filters out points which don’t fall into the specified range of Z values. We can use these outliers to limit the range of values (748-970 m seems to be a safe bet). Use the zrange parameter to filter out any possible outliers (don’t be confused with the intensity_range parameter) when computing a new DSM:

r.in.lidar input=input=297285v.las output=binned_dsm_limited method=max -o zrange=748,970

Interpolation

Now we will interpolate a digital surface model (DSM) and for that we can increase the resolution to obtain as much detail as possible (we could use 1, i.e. res=1 foot, for high detail or 3, i.e. res=3 feet -a, to increase the computational speed). I will use 3 feet, but feel free to play with the resolution:

g.region raster=Arb_297285 res=3 -ap

Before interpolating, let’s confirm that the spatial distribution of the points allows us to safely interpolate. We need to use the same filters as we will use for the DSM itself:

r.in.lidar input=input=297285v.las output=count_dsm method=n return_filter=first -o zrange=748,970

Here we evaluate if there are possible outliers:

r.report map=count_dsm units=h,c,p

Then check the spatial distribution. For that we will need to use histogram equalized color table (legend can be limited just to a specific range of values: d.legend raster=count_interpolation range=0,5).

r.colors map=count_dsm color=viridis -e

First we will import the LAS file as vector points. We will only use the first return points, and limit the import vertically to avoid using outliers found in the previous step. Before running it, uncheck Add created map(s) into layer tree in the v.in.lidar dialog if you are using GUI.

v.in.lidar -t -b -o input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=first_returns zrange=748,970

Then we interpolate (if you have problems with the calculation you may need to set the extent using g.region raster=dsm):

v.surf.rst input=first_returns elevation=dsm tension=25 smooth=1 npmin=80

Terrain analysis

LiDAR can also be used to identify different terrain features. In this exercise we will identify different physical features of the landscape.

If you haven’t already, let’s set the computational region based on an existing raster map:

g.region raster=dsm -p

We will now import points based on the ground points only (the -t flag disables creation of attribute table and the -b flag disables building of topology; uncheck Add created map(s) into layer tree):

v.in.lidar -bt input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=points_ground class_filter=2 -o zrange=748,970

This figure indicates the different classes of the points in the LiDAR point cloud:

Classification value Meaning
0 Never classified
1 Unassigned
2 Ground
3 Low Vegetation
4 Medium Vegetation
5 High Vegetation
6 Building
7 Low Point
8 Reserved*
9 Water
10 Rail
11 Road Surface
12 Reserved*
13 Wire - Guard (Shield)
14 Wire - Conductor (Phase)
15 Transmission Tower
16 Wire-Structure Connector (Insulator)
17 Bridge Deck
18 High Noise
19-63 Reserved
64-255 User Definable

Now we will interpolate the ground surface using regularized spline with tension (implemented in v.surf.rst) and at the same time we also derive slope, aspect and curvatures (the following code is one long line):

v.surf.rst input=points_ground tension=25 smooth=1 npmin=100 elevation=terrain slope=slope aspect=aspect pcurvature=profile_curvature tcurvature=tangential_curvature mcurvatur=mean_curvature

When we examine the results, especially the curvatures show a pattern which may be caused by some problems with the point cloud collection. We decrease the tension which will cause the surface to consider fewer points and increase the smoothing between points. Since the raster map called range already exists, use the overwrite flag, i.e. –overwrite in the command line or checkbox in GUI to replace the existing raster by the new one. We use the –overwrite flag shortened to just –o.

v.surf.rst input=points_ground tension=20 smooth=5 npmin=100 elevation=terrain slope=slope aspect=aspect pcurvature=profile_curvature tcurvature=tangential_curvature mcurvatur=mean_curvature --o

When we are satisfied with the result(feel free to play with the tension and smooth variables), we can compute shaded relief:

r.relief input=terrain output=relief

Now combine the shaded relief raster and the elevation raster to create a shaded relief map. Shaded relief, or hill-shading, shows the shape of the terrain in a realistic fashion by depicting how the three-dimensional surface would be illuminated from a point light source. We can calculate relief in several ways. In the GUI, we produce shaded-relief by combining the layers and changing opacity of layer or the other. We can also do this programmatically using the r.shade module, which combines the two rasters and creates a new one. We can also do this on the fly without creating any raster using the d.shade module. The module can be used from the GUI through the toolbar or using the Console tab:

d.shade shade=relief color=terrain

Alternatively, we can calculate the skyview (how it looks from an airplane) using r.skyview:

r.skyview input=terrain output=skyview ndir=8 colorized_output=terrain_skyview

Let’s combine the terrain and skyview on the fly:

d.shade shade=skyview color=terrain

We can also create a shaded relief map that highlights terrain features which wouldn’t be visible using standard shading technique. The module show relief from various directions and combines them into RGB composition:

r.shaded.pca input=terrain output=pca_shade

Local relief model. Local relief models enhance the visibility of small-scale surface features by removing large-scale landforms from the DEM:

r.local.relief input=terrain output=lrm shaded_output=shaded_lrm

The r.terrain.texture module can be used to calculate different terrain characteristics. See the the help file here. For example, we can calculate pit density:

elevation=terrain thres=0 pitdensity=pit_density

Finally, we can use the r.geomorphon module, which automatically detection of landforms (using 50 m search window):

-m elevation=terrain forms=landforms search=50

Vegetation analysis

We can also isolate vegetation using LiDAR. To do this we will use the zrange parameter to filters out points which don’t fall into the specified range of Z values. Although for many vegetation-related application, a coarse resolution is necessary/appropriate due to the number of points needed to interpolate the map, we will use 1 meter:

g.region raster=dsm res=3

Let’s start by exploring the range of all the points by estimating heights for each cell. This allows us to identify tall objects infered from the difference between points at a single location:

r.in.lidar input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=range method=range -o zrange=748,970

To isolate and visualize tall objects, we can use r.mapcalc. This is a raster calculator in Grass. All expressions are noted using quotation marks. Here we are using a logical operator if to identify in the range map locations that are greater than 2 meters (6 feet). The cells or raster with a range greater than 2 meters are given 1, and those that are less than 2 meters are given as null:

r.mapcalc "vegetation_by_range = if(range > 6, 1, null())"

We can change the color table for a raster map from layer context menu using Set color table interactively.

Now let’s calculate the tallest objects in relation to the actual terrain. We do this by using the max points and using the base_raster we already calculated. The function assumes that the base map is zero and we set the range to 100 meters as it is unlikely that there are such tall objects

r.in.lidar -d input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=height_above_ground method=max base_raster=terrain zrange=0,100 -o

Now let’s redo the calculation to identify locations with objects that are higher than 2 m. If we wanted to preserve the elevation in the output, we could use if(height_above_ground > 6, height_above_ground, null()) in the expression, but we want just the areas, so we use constant (1).

r.mapcalc "above_2m = if(height_above_ground > 6, 1, null())"

Often there are speckles where there is no data due to issues with the point cloud interpolation. To reduce these speckles, we can use the r.grow, which assign data around existing non-null cells/rasters:

r.grow input=above_2m output=vegetation_grow

This is likely a good approximation of locations with vegetation and the occasional building (which we can filter out later). We calculate the clump of individual cells into patches using r.clump. This helps identify contiguous area of forests and filter out small insignificant vegetation:

r.clump input=vegetation_grow output=vegetation_clump

Some of the patches are very small. Using r.area we remove all patches smaller than the given threshold:

r.area input=vegetation_clump output=vegetation_by_height lesser=100

Now convert these areas to vector:

r.to.vect -s input=vegetation_by_height output=vegetation_by_height type=area

So far we were using elevation of the points, now we will use intensity. Intensity is a measure, collected for every point, of the return strength of the laser pulse that generated the point. It is based, on the reflectivity of the object struck by the laser pulse intensity is used by r.in.lidar when -j (or -i) flag is provided:

r.in.lidar input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=intensity zrange=748,970 -j -o

Given this high resolution, there are some cells without any points, so we use r.fill.gaps to fill these cells (r.fill.gaps also smooths the raster as part of the gap-filling process).

r.fill.gaps input=intensity output=intensity_filled uncertainty=uncertainty distance=3 mode=wmean power=2.0 cells=8

Let’s apply a grey scale color table as it is more appropriate for intensity:

r.colors map=intensity_filled color=grey

There are some areas with very high intensity. We can equalize the color table using the histogram to make it easier to visually explore the data

r.colors -e map=intensity_filled color=grey

Let’s use r.geomorphon again, but now with DSM to identify specific features that can be inferred from the intensity layer. The following settings shows building roof shapes:

r.geomorphon elevation=dsm forms=dsm_forms search=7 skip=4

With different settings, especially a higher flatness threshold (flat parameter), we can isolate forested areas. Using the geomorphon function vegetation appears as footslopes, slopes, and shoulders. Individual trees are represented as shoulders.

r.geomorphon elevation=dsm forms=dsm_forms search=12 skip=8 flat=10 --o

Lowering the skip parameter decreases generalization and brings in summits which often represents tree tops:

r.geomorphon elevation=dsm forms=dsm_forms search=12 skip=2 flat=10 --o

Now we extract summits (trees) using r.mapcalc raster algebra expression if(dsm_forms==2, 1, null()).

r.mapcalc "trees = if(dsm_forms==2, 1, null())"

How would you get a vectors of tree locations?

3D geovisualization of rasters

We can explore our study area in the 3D view:

alt text alt text

We might also want to mask the data to a specific location using r.mask(this is clip in many gis platforms). You can download .las file for Ann Arbor here() for doing your own terrain models. Ensure that you choose the correct computational surface.

Optional

Basic introduction to Python interface

You can execute Python code in the Simple Python Editor (accessible from the toolbar or the Python tab in the Layer Manager). Another option is to use your favorite text editor and then run the script in GRASS GIS using the main menu File → Launch script.

We can use the Simple Python Editor to run commands. When you opening the Editor, you find a short code snippet that imports the GRASS GIS Python Scripting Library:

import grass.script as gscript

You can run Grass function using the gscript.run_command. For example, if we want to call the g.region function we run this code:

gscript.run_command('g.region', flags='p')

In this example, we set the computational extent and resolution to the raster layer ortho which would be done using g.region raster=ortho in the command line. To use the run_command to set the computational region, replace the previous g.region command with the following line: